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FROM DYNAMICS TO LINKS: A SPARSE RECONSTRUCTION OF THE 
TOPOLOGY OF A NEURAL NETWORK 


GIACOMO ALETTI, DAVIDE LONARDONI, GIOVANNI NALDI, AND THIERRY NIEUS 

ABSTRACT. One major challenge in neuroscience is the identification of interrelations be¬ 
tween signals reflecting neural activity and how information processing occurs in the neural 
circuits. At the cellular and molecular level, mechanisms of signal transduction have been 
studied intensively and a better knowledge and understanding of some basic processes of 
information handling by neurons has been achieved. In contrast, little is known about the 
organization and function of complex neuronal networks. Experimental methods are now 
available to simultaneously monitor electrical activity of a large number of neurons in real 
time. Then, the qualitative and quantitative analysis of the spiking activity of individual 
neurons is a very valuable tool for the study of the dynamics and architecture of the neural 
networks. Such activity is not due to the sole intrinsic properties of the individual neural 
cells but it is mostly consequence of the direct influence of other neurons. The deduc¬ 
tion of the effective connectivity between neurons, whose experimental spike trains are 
observed, is of crucial importance in neuroscience: first for the correct interpretation of the 
electro-physiological activity of the involved neurons and neural networks, and, for cor¬ 
rectly relating the electrophysiological activity to the functional tasks accomplished by the 
network. In this work we propose a novel method for the identification of connectivity of 
neural networks using recorded voltages. Our approach is based on the assumption that the 
network has a topology with sparse connections. After a brief description of our method 
we will report the performances and compare it to the cross-correlation computed on the 
spike trains, that represents a gold standard method in the field. 


1. Introduction 

Along the latests years there have been enormous progresses in the Neuroscience field 
that have revolutionized the way we envisage the brain functions. In this regard, the tech¬ 
nological advancements have been fundamental to improve the recording capability from 
brain areas and neural populations El. Nowadays, multi-site recordings can be achieved 
from thousands of channels (sites) with a good spatial (at the cellular level) and temporal 
resolution (less than one millisecond for the action potential) yielding a good description of 
the underlying network dynamics. Given that the brain operates on a single trial basis such 
recordings are instrumental to understand the neural code a. As a first step, multi-site 
recordings allow to quantify the information flow in the network. The anatomical wiring 
(i.e. Structural Connectivity, SC) clearly plays a fundamental role to understand how cells 
comunicate among them but it is often not well known neither it can by itself explain the 
overall network activity. Multi-site recordings can be used to infer statistical dependencies 
(i.e. Functional Connections, FC) among the recorded units and to track the information 
flow in the network £3). On the other hand the Effective Connectivity (EC) denotes the 
directed causal relationship between the recorded sites. The EC is typically estimated by 
stimulating one cell and studying the effects on the connected elements. Alternatively the 
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EC can also be studied establishing a causal mathematical model between the recorded 
units data. 

Importantly, multi-site recordings raise some limitations that need to be evaluated care¬ 
fully before any further analysis. First, the experimental sessions are often severly limited 
in time. Second, the high dimensional data sets involve a set of numerical and mathemati¬ 
cal problems that would be hard to face even with long enough recording sessions. These 
issues are common to different fields and have been coined as ’curse of dimensionality’. To 
overcome these issues, two approaches are typically foreseen. A first solution consists into 
a reductionist approach that projects the data on a lower dimenstional space that can better 
elucidate the underlying processing. Another possibility consists into fitting the observed 
data to a low-dimensional model that captures the salient properties of the dynamics EG). 
Here we introduce a model of effective connectivity that gets rid of the dimensionality 
problem by introducing a natural constrain of almost all biological networks: the sparse¬ 
ness among the connected units (see e.g. fT|). Other approaches, such as the multivariate 
autoregressive model ESI, allow for the possibility of a fully connected network in which 
every node may influence all the other nodes. However this is somewhat unrealistic for 
biological networks (i.e. each neuron is directly influenced by only a small subset of neu¬ 
rons) and it also leads to practical challenges. In fact, most of the connectivity tools used 
to understand the communication among neuronal populations are based on linear models 
although it is widely recognized that the interactions (i.e. synaptic currents) among neural 
cells are nonlinear. Then, fully connected networks involve models that can easily be over¬ 
parameterized and hard to solve. Therefore current methods are not well suited to robustly 
infer on the network connectivity from the time series in different contexts. Among these 
approaches, the Granger causality HE) is probably the most prominent and most widely 
used concept. This concept of causality does not rely on the specification of a scientific 
model and thus is particularly suited for empirical investigations of cause-effect relation¬ 
ships. On the other hand, if important relevant variables are not included in the analysis, 
the Granger causality can lead to so-called spurious causalities . In order to capture non¬ 
linear interactions between even short and noisy time series, we consider an event-based 
model. Then, we involve the physiological basis of the signal, which is likely non-linear. 

In Section [2] we introduce a general setting for the problem. In Section [3] we describe 
our experimental settings and data preprocessing. The results are reported in Section [4] 
while final remarks are postponed in the last section. 

2. Problem Formulation 

2.1. Network representation of the problem. Let us consider a graph G = (V, E), 

where V corresponds to a set of N neurons connected with a directed graph given by 
the edges contained in E of ordered pairs of vertices contained in V x V. On each node 
i £ V, two stochastic processes are observed: 

• an external process {Xi(t), 0 < t < T} (spike); 

• an internal process 0 < t < T}. 

We assume a Local Markov Property (LMP) : the relevant information on each process 
Yi(t),i = 1,.... A is “contained” on the external activity of the neurons j connected to 
it, and in a time interval S. For example, in the network given in Figure [T] the law of the 
4th neuron depends on the external activity of the 1st, 2nd, and 5th neuron together with 
its own internal activity. 

We explain now the LMP in terms of the er-algebras of events observable by the pro¬ 
cesses. We denote by T f - the observable events of the “past” of the external processes 
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until time t, and it will be defined as the a- algebra of events generated by the processes 
{Xi(s),s < t,i e V}: T t - = Viev cr(Xi(s),s < t). The u-algebra contains 
the observable events of the “close past” of the in-neighborhoods of the i-th neuron: 
Xfl = V, • (j i)&E a (Xj{s),t — 5 < s < t). Moreover, we will denote by Q t + the 

observable events of the “future” of the internal processes: = a(Yi(u),t < u). 

With this notation, LMP says that the future of the i-th internal process and past of all 
the external processes are conditionally independent given the close past of the i-th neuron: 

(LMP) P{G\F t -) = P(G\J$>_), VGeSg.Vi. 

(i) 

In other words, J- & gives all the relevant information on the future of the internal activity 
of the i-th neuron. 



/0 1 0 1 0 \ 

0 0 0 1 0 

1 0 0 0 1 

0 0 10 0 
\1 0 0 1 0 / 


FIGURE 1. Example of a nework of 5 neurons and its adjacency matrix. 


2.2. Event-based discrete time model. We use a particular model of the class of pro¬ 
cesses satisfying ( |LMP| i. We sample the processes at discrete times t m = m5 , and we 
assume an event-based model : both the internal and external processes are simple point 
processes, and we assume that at most one event may occur for each process in any time 
interval. We note that the definition of event depends on the knowledge of the investigator 
and on the available data. For the purpose of this paper, the definition of events will be 


introduced in the Section 3.2 Given an event-based model, for any m = 1,..., M, we 
observe 


Vi = 


+1 if there has been an event for Yj (s) duting (t m - 1 , t m ]; 
0 otherwise; 


and 


+1 if there has been an event for Xi(s) duting t m \; 

0 otherwise; 


The model can be completely characterized by defining the family of conditional probabil¬ 
ities: for any i, m. 


7t(*i m) = P(y z m+1 = 1 \x l j,l< m,j eV) = P(y™ +1 = l\x™,(j,i) £ E) 

the last equality being a consequence of ( |LMP| ). 

To model the association among the nodes, we assume in this paper a time-homogeneous 
linear logit model: 


(1) = (1 + exp(—/3 0 — ^ Pj,iX™)) 1 
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Given a set {cui,m} of positive real numbers, the negative weighted log-likelihood function 
for our processes reads: 

M—l 

(2) £ m (J3) = - Y 5Z Wi - m (^ m+ll °g 7r (* I TO ) + (! - y» T " +1 ) 1 °g( 1 - 7r(i, m))) 

m= 1 

Note that the significant parameters {/3j i, (j, i) £ E} of our model depend on the topology 
of the network. The choice of the weights is done here to balance the number of y" L = 0 
and y™ = 1. A simple choice might be uj^ m = i v’j '^VT = 0 and w i)m = ;( 1_ 2/j) 

if»r = i- 

3. Experimentimental settings 

3.1. Simulations . We simulate the potential of the j -th neuron Vj (t) by the Hodgkin&Huxley 
equations of the cerebellar granule cell (GrC, Q). Heterogeinity among the GrCs is in¬ 
troduced by randomizing the resting potentials with a current of 2 ± 0.2 pA. The network 
consists of 20 neurons with a connectivity probability of 0.3 and all connections are direc¬ 
tional (i.e. chemical synapses). Since self-connections (i.e. ’autapses’) are not so frequent 
in the brain, they were not included in our network. Based on these numbers, the simulated 
network comprises on average 114 connections (i.e. 20 • 19 • 0.3). The chemical synapses 
are modeled by the fast AMPA excitatory currents IfTOl . In addition, to mimic a biologi¬ 
cally realistic noisy regime, Poisson distributed AMPA currents (of frequency 0.2 Hz) are 
injected to all GrCs. The cellular and synaptic models are described in the literature Sosa 
and the parameter settings and their changes respect to the literature are reported in the 
Table [I] Since here we are interested in highlighting the potential impact of our approach 


synaptic input 

parameter 

unit 

value 

value in literature [11011 

AMPA 

gmax 

pS 

800 

1200 


n 

ms l mM 1 

5.4 



7*2 

ms -1 

0.84 



r& 

1 

0 

1.12 

NMDA 

not included here 




synaptic noise 

gmax 

pS 

500 



n 

ms~ l mM ~ 1 

5.4 



T2 

ms~ l 

0.1 



re 

1 

0 



Table 1. Values for the parameter settings for the dynamics of AMPA, 
NMDA and synaptic currents 


we discarded the contribution of other synaptic conductances that will be included in a next 
work. 

We simulate five seconds of activity of such a network. A snapshot of the activity of a 
neuron in the network is shown in Figure[2] The upswings of the potential are mainly deter¬ 
mined by the inputs from the other cells (the time occurrence of these inputs is highlighted 
by the red, blue and green dashed lines). Interestingly, the spiking activity is not strictly 
determined by the input itself. In fact, the spike doublet (D1,D2), the isolated spikes (IS) 
and the spike from excitation (EXC) are not apparently determined by the input. Spike 
D1 arises from a depolarization that is not directly caused by the input and spike D2 can 
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either by determined by high membrane excitation as well as by noise. The spike IS is 
most probably caused only by noise and spike EXC is again due to membrane excitation 
or simply by synaptic noise. 



Figure 2. Voltage trace of a neuron in the network. The voltage (black 
trace) is overlayed with the time stamps of the inputs (colored dash lines) 
to the cell. The voltage changes are clearly correlated with the input. 
However spikes Dl, D2, IS, EXC are not apparently determined solely 
by the input. Note that the blue, red and green colors correspond to three 
different cells. 


3.2. Data preprocessing. The process {xj 1 , m = 1, is the discretization (in 

time) of the spike activity of the j-th neuron, so that x™ = 1 if there has been a spike 
during the time interval i m ]. 

The process {yj 1 , m = 1,..., M} is a reaction process, nonlinear filter of the potential 
activity Vj(t). Here, y 7 - 1 = 1 if there has been an event for Vj(t) during the time interval 
t m ], where an event at time t depends on three factors (see Figure [3ji: 

(i) the right derivative V- (t + ) must be greater than a positive threshold (increasing of 
potential after excitation); 

(ii) the increasing of derivatives V- (£+) — V-(t~) must be greater than a positive 
threshold (convex effect of potential due to excitation); 

(iii) the left derivative VJ(t~ ) must be greater than a negative threshold (in connection 
with other conditions, this avoid an event caused by the recover of the resting 
potential during the hyperpolarization phase; alone, this identifies those events). 

More precisely, when a sequence of times of length n > 1 satisfies the requested con¬ 
ditions, an event is detected. The first time of this sequence is the time of the event. 

3.3. Topology estimator. Given the data of our processes, we are interested in recon¬ 
structing the topology of our network. We want to give an estimator V of the set V of the 
edges. In our contest, the edge (j. i) exists when the process y-" is directly caused by x™, 
in the sense of ( |LV1 F[ >. In other words, when f3jy is different from zero in Q- Therefore, 
we adopt the following strategy: 
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Figure 3 . Example of identification of events based on the voltage ac¬ 
tivity. Events that are determined by the links of the network are marked 
with * (i.e. true positive), while o shows false identification (i.e. false 
positive) (due, for example, to external noise). Top left: identification of 
events with all the three factors of data preprocessing. Top right: identifi¬ 
cation of events with only (i) factor. Bottom left: identification of events 
with only (ii) factor. Bottom right: identification of events with only (iii) 
factor. 


• with a penalization technique, we find a sparse estimator (3 of /3; 

• we say that (j, i) £ V if Pj,i is different from zero. 

In this paper, we adopt a t \-penalization on the regression coefficients i,j £ E}. 
Given a positive penalization parameter A, let 

HP, A) = M/3) + A ^ \fo\, 

i,j£E 


and define the Lasso estimator /3(A) = argmin^ L(/3, A). As stated above, we define 
V{\) as 

(3) {j,i) e V <=>• /3j,j(A) > 0. 


Here we also compare the performances of the novel methodology with the standard 
cross-correlation that is widely used in the multi electrode array field 0(9: T31. The cross¬ 
correlation functions among the discrete spike trains are defined as: 


(ST t • STj(r)) 


(4) 


CC tJ (r) 
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where ST,;, 57} are the binned spike trains (bin size 1 ms) and Ni,Nj the corresponding 
number of spikes. The strength of a connection, between the nodes i and j, is then given 
by the peak of the quantity CCij given in <[4]». The network topology can then be inferred 
by retaining the strongest and most significants cross-correlation peaks that are overcome 
a selected threshold. 


4. Experimental Results 


The validation of the proposed methodology is afforded on simulations presented in the 
The performances of the methodologies are quantified with the well known 


Section 3.1 


receiver operating characteristic (ROC) curves and with the recently introduced positive 
precision curve (PPC, ©). The graphs reported in Figure [4] are obtained by varying the 
penalization parameter A, for the Lasso methodology, and the cross-correlation threshold, 
for the respective correlation methodology. 

Let us remind that PPC is defined as: 


(5) 


PPC 


TP-FP 
TP+ FP 


and represents the proportion of the correctly (true positive, TP) versus the incorrectly 
(false positive, FP) inferred links. The plot on the right in the Figure [4] reports the PPC 
index with respect to the number of links included in the analysis. Interestingly, the PPC 
curve of the Lasso stays at its maximum up to 30% of the links, that corresponds to the 
real links of the network (the connectivity probability is 0.3). Clearly, beyond 30% of the 
included links, the PPC has a negative power decay. We point out that the cross-correlation 
does not achieve to infer the topology at any level of the cross-correlation threshold. In¬ 
terestingly, when the Lasso events are based only on spiking information (condition (iii) 
alone) it performs as bad as the cross-correlation. The Lasso (Figure^ performs very well 
and reaches the maximum value (=1) in both the ROC and PPC curves for a large range 
of the penalization parameter. Moreover, the worst performances are achieved when no 
constraint is imposed on the sparsity (e.g. A = 0) of the inferred network. Finally, the 
intuition is preserved: the first connections introduced in the Lasso estimate, are almost all 
true (PPC ~ 1). Another interesting feature of the newly introduced Lasso methodology 
consists into its capability of inferring the birectional links. The results of Figure [4] are 
based on a network with 121 connections out of which 38 were bidirectional (31.4% out 
of the total). This is an another advantage over the cross-correlation method, that only 
determines unidirectional connections. 


5. Conclusions and Future Works 

Understanding how interactions between brain structures support the performance of 
specific cognitive tasks or perceptual processes is a prominent goal in neuroscience. The 
effects that one part of the nervous system has on another are typically examined by stimu¬ 
lating or lesioning the first part and investigating the outcome in the second. For example, 
in peripheral and spinal pathways, the interventional techniques of stimulation and abla¬ 
tion have proven to be powerful methods for inferring causal influences from one neuron 
or neuronal population to another. For the study of causal relations within the brain (func¬ 
tional and effective connectivity ), the utility of the interventional techniques is diminished 
by the high levels of convergence and divergence in brain pathways. 

In this work we have developed an event based approach for inferring networks of causal 
relationships in a neuronal population. Specifically, we suppose that we are able to observe 
the dynamical behaviors of individual components of a neuronal networks and that few of 
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FIGURE 4. ROC (on the left) and PPC (on the right) curves for the topol¬ 
ogy estimator with the complete data preprocessing and with each 
of the filtering conditions (i), (ii) and (iii) given in Section 3.2 The first 
conditions (i) and (ii) of data preprocessing may be used alone without 
affecting the main result, while spiking activity alone is not predictive 


the components may be causally influencing each other. The variables could be invasive 
electrode recordings, intracranial EEG, or non-invasive EEG, MEG or fMRI time series 
from different parts of the brain. In order to introduce our method we have considered a 
simulated cerebellar granule cell network capturing nonlinear interactions between even 
short and noisy time series. 

The results we got are quite promising from many point of views. First, despite the 
algorithm was applied only on short simulations (5 seconds) it achieved to filter out noisy 
from causal responses yielding a reliable estimate of the underlying connections in the net¬ 
work. Second, the approach is quite general and the conditions (i),(ii),(iii) can be further 
adapted to different types of electro-physiological signals. Third, the Lasso is also quite 
robust respect to bidirectional connections in the network. This is of fundamental impor¬ 
tance since bidirectional network motifs are quite abundant in the brain ED. The proposed 
Lasso methodology assumes the knowledge of the voltage traces. From an experimental 
point of view such a detailed information can be achieved with patch-clamp experiments 
but this technique is not designed to perform simultaneous recordings from populations 
of neurons. Interestingly, nowadays there have been huge improvements in the field of 
the genetically encoded voltage indicators that allow to track sub-threshold activities |7), a 
key ingredient of the Lasso algorithm. These progresses will likely provide in a near future 
multi-site recordings of a population of neurons. 

As a next step, we will then test the robustness of the methodology including differ¬ 
ent noise sources (i.e. membrane noise), inhibitory connections and verify its robustness 
respect to bigger networks. 


References 

[1] A. Bolstad, B. Van Veen, and R. Nowak. Causal network inference via group sparse regularization. IEEE 
Transactions on Signal Processing , 59(6):2628-2641, 2011. 

[2] S. L. Bressler and A. K. Seth. Wienergranger causality: A well established methodology. Neuroimage , 

58(2):323329,2011. ^ 

[3] E. Bullmore and O. Sporns. Complex brain networks: Graph theoretical analysis of structural and functional 
systems. Nature Reviews Neuroscience , 10(3): 186-198, 2009. 

[4] M. M. Churchland, B. M. Yu, M. Sahani, and K. V. Shenoy. Techniques for extracting single-trial activity 
patterns from large-scale neural recordings. Curr Opin Neurobiol , 17(5):609-618, Oct 2007. 

















FD2L: SPARSE RECONSTRUCTION OF THE TOPOLOGY 


9 


[5] E. D’Angelo, T. Nieus, A. Maffei, S. Armano, P. Rossi, V. Taglietti, A. Fontana, and G. Naldi. Theta- 
frequency bursting and resonance in cerebellar granule cells: Experimental evidence and modeling of a 
slow k+-dependent mechanism. Journal of Neuroscience , 21(3):759-770, 2001. 

[6] M. Garofalo, T. Nieus, P. Massobrio, and S. Martinoia. Evaluation of the performance of information theory- 
based methods and cross-correlation to estimate the functional connectivity in cortical networks. PLoS One , 
4(8):e6482, 2009. 

[7] Y. Gong, C. Huang, J. Z. Li, B. F. Grewe, Y. Zhang, S. Eismann, and M. J. Schnitzer. High-speed record¬ 
ing of neural spikes in awake mice and flies with a fluorescent voltage sensor. Science (New York, N.Y.) , 
350:1361-1366, Dec 2015. 

[8] C. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica , 
37:424-438, 1969. 

[9] A. Maccione, M. Garofalo, T. Nieus, M. Tedesco, L. Berdondini, and S. Martinoia. Multiscale functional 
connectivity estimation on low-density neuronal cultures recorded by high-density cmos micro electrode 
arrays. J Neurosci Methods , 207(2): 161—171, Jun 2012. 

[10] T. Nieus, E. Sola, J. Mapelli, E. Saftenku, P. Rossi, and E. D’Angelo. Ltp regulates burst initiation and 
frequency at mossy fiber-granule cell synapses of rat cerebellum: experimental observations and theoretical 
predictions. J Neurophysiol , 95(2):686-699, Feb 2006. 

[11] S. Song, P. J. Sjostrom, M. Reigl, S. Nelson, and D. B. Chklovskii. Highly nonrandom features of synaptic 
connectivity in local cortical circuits. PLoS Biol , 3(3):e68, Mar 2005. 

[12] I. H. Stevenson and K. P. Kording. How advances in neural recording affect data analysis. Nat Neurosci , 
14(2): 139-142, Feb 2011. 

[13] S. Ullo, T. R. Nieus, D. Sona, A. Maccione, L. Berdondini, and V. Murino. Functional connectivity estima¬ 
tion over large networks at cellular resolution based on electrophysiological recordings and structural prior. 
Front Neuroanat , 8:137, 2014. 

[14] M. Winterhalder, B. Schelter, W. Hesse, K. Schwab, L. Leistritz, D. Klan, R. Bauer, J. Timmer, and H. Witte. 
Comparisson of linear signal processing techniques to infer directed interactions in multivariate neural sys¬ 
tems. SignalProcess^ 85(11):21372160, 2005. 

ADAMSS Center, Universita degli Studi di Milano, 20131 Milano, Italy 

E-mail address : giacomo . aletti@unimi . it 

Neuroscience Brain Technology, Istituto Italiano di Tecnologia, via Morego 30, 16163 

Genova, Italy 

ADAMSS Center, Universita degli Studi di Milano, 20131 Milano, Italy 

Neuroscience Brain Technology, Istituto Italiano di Tecnologia, via Morego 30, 16163 

Genova, Italy 



